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Numerical integration of the differential equations of light propagation in the 
Schwarzschild metric shows that in some extreme situations relevant for practical 
observations (e.g. for Gaia) the well-known standard post-Newtonian formula for 
the boundary problem has an error up to 16 /uas. The aim of this note is to identify 
the reason for this error and to derive an extended formula accurate at the level of 
1 /ias as needed e.g. for Gaia. 

The analytical parametrized post-post-Newtonian solution for light propagation 
derived by Klioner & Zschocke [l(| gives the solution for the boundary problem with 
all analytical terms of order 0{c~'^) taken into account. Giving an analytical upper 
estimates of each term we investigate which post-post-Newtonian terms may play a 
role for an observer in the solar system at the level of 1 /ias. We conclude that only 
one post-post-Newtonian term remains important for this numerical accuracy and 
derive a simplified analytical solution for the boundary problem for light propagation 
containing all the terms that are indeed relevant at the level of 1 /ias. The derived 
analytical solution has been verified using the results of a high-accuracy numerical 
integration of differential equations of light propagation and found to be correct at 
the level well below 1 /tas for arbitrary observer situated within the solar system. 



2 



Contents 

I. Introduction 

II. Schwarzschild metric and null geodesies in harmonic coordinates 

A. Metric tensor 

B. Christoffel symbols 

C. Isotropic condition 

D. Equation of isotropic geodesies 

III. Numerical Integration of the equations of light propagation 

IV. Standard post-Newtonian approach 

A. Equations of post-Newtonian approach 

B. Comparison between the post-Newtonian approximation and numerical solution [8 

V. Post-post-Newtonian solution of boundary problem ^ 

A. Analytical estimates of the individual terms in cr 10 

B. Analytical estimates of the individual terms in transformation from fc to cr 11 

C. Analytical estimates of individual terms in transformation from cr to n 13 

D. Analytical estimates of individual terms in transformation from k to n [14 



VI. Transformation from fc to n for stars and quasars 

A. Transformation from fc to cr 

B. Transformation from cr to n 



ii 

12 

1^ 



VII. Summary and concluding remarks |20 



Acknowledgments |21 



References 21 



A. Estimates of terms in the Shapiro delay |22 



B. Estimates of terms in the transformation between cr and fc 22 



C. Estimates of terms in the transformation between n and cr 24 



D. Estimates of terms in the transformation between n and fc 24 



E. Estimates of terms in the transformation between n and cr for stars and 
quasars 



26 



3 



I. INTRODUCTION 

It is well known that adequate relativistic modelling is indispensable for the success of 
microarcsecond space astrometry. One of the most important relativistic effects for astro- 
metric observations in the solar system is the gravitational light deflection. The largest 
contribution in the light deflection comes from the spherically symmetric (Schwarzschild) 
parts of the gravitational fields of each solar system body. Although the planned astrometric 
satellites Gaia, SIM, etc. will not observe very close to the Sun, they can observe very close 
to the giant planets also producing significant light deflection. This poses the problem of 
modelling this light deflection with a numerical accuracy of better than 1 fias. 

The exact differential equation of motion for a light ray in the Schwarzschild field can be 
solved numerically as well as analytically. However, the exact analytical solution is given 
in terms of elliptic integrals, implying numerical efforts comparable with direct numerical 
integration, so that approximate analytical solutions are usually used. In fact, the stan- 
dard parametrized post-Newtonian (PPN) solution is sufficient in many cases and has been 
widely applied. So far, there was no doubt that the post-Newtonian order of approximation 
is sufficient for astrometric missions even up to microarcsecond level of accuracy, besides as- 
trometric observations close to the edge of the Sun. However, a direct comparison reveals a 
deviation between the standard post-Newtonian approach and the exact numerical solution 
of the geodetic equations. In particular, we have found a difference of up to 16 /las in light 
deflection for solar system objects observed close to giant planets. This error has triggered 
detailed numerical and analytical investigations of the problem. 

Usually, in the framework of general relativity or the PPN formalism analytical orders of 
smallness of various terms are considered. Here the role of small parameter is played by 
where c is the light velocity. Standard post-Newtonian and post-post-Newtonian solutions 
are derived by retaining terms of relevant analytical orders of magnitude. On the other hand, 
for practical calculations only numerical magnitudes of various terms are relevant. In this 
note we attempt to close this gap and combine the analytical post-post-Newtonian solution 
derived in Klioner & Zschocke p| with estimates of numerical magnitudes of various terms. 
In this way we will derive a compact analytical solution for the boundary problem for light 
propagation where all terms are indeed relevant at the level of 1 /xas. The derived ana- 
lytical solution is then verified using high-accuracy numerical integration of the differential 
equations of light propagation and found to be correct at the level well below 1 fias. 

We use fairly standard notations: 

• G is the Newtonian constant of gravitation; 

• c is the velocity of light; 

• (3 and 7 are the parameters of the Parametrized Post-Newtonian (PPN) formalism 
which characterize possible deviation of the physical reality from general relativity 
theory (/3 = 7 = 1 in general relativity); 

• lower case Latin indices i, j, . . .take values 1, 2, 3; 

• lower case Greek indices /z, z/, . . .take values 0, 1, 2, 3; 

• repeated indices imply the Einstein's summation irrespective of their positions (e.g. 

a' V = a^h^ + 6^ + ^3 ^3 = + + a'^ b'^ + b^); 
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• a dot over any quantity designates the total derivative with respect to the coordinate 

time of the corresponding reference system: e.g. a = — ; 

at 

• the 3-dimensional coordinate quantities ( "3- vectors" ) referred to the spatial axes of 
the corresponding reference system are set in boldface: a = a*; 

• the absolute value (Euclidean norm) of a "3- vector" a is denoted as \a\ or, simply, a 
and can be computed as a = \a\ = (a^ + + a^)^/^; 

• the scalar product of any two "3-vectors" a and b with respect to the Euclidean metric 
6ij is denoted by a ■ 6 and can be computed as a ■ b = 6ij a* b' = b^; 

• the vector product of any two "3-vectors" a and b is designated by a x 6 and can 
be computed as (a x b)* = Eijta^ where Eijk = {i — — k){k — i)/2 is the fully 
antisymmetric Levi-Civita symbol; 

• for any two vectors a and b the angle between them is designated as 5(a, b). 

The paper is organized as follows. In Section[TTlwe present the exact differential equations. 
High-accuracy numerical integration of these equations is discussion in Section HTTl In Section 
IIVI we discuss the standard post-Newtonian approximation and demonstrate the problem 
with the standard post-Newtonian solution by direct comparison between numerical results 
and the PPN solution. In Section |Vl the formulas for the boundary problem in post-post- 
Newtonian approximation are considered. A detailed estimation of all relevant terms is given, 
and simplified expressions are derived. We demonstrate by explicit numerical examples the 
applicability of this analytical approach for the GAIA astrometric mission. In Section IVII 
we consider the important case of objects situated infinitely far from the observer as a limit 
of the boundary problem. The results are summarized in Section IVIII In the Appendices 
detailed derivations for a number of analytical formulas are given. 

II. SCHWARZSCHILD METRIC AND NULL GEODESICS IN HARMONIC 

COORDINATES 

For the reasons given above we need a tool to calculate the real numerical accuracy 
of some analytical formulas for the light propagation. To this end, we consider the exact 
Schwarzschild metric and its null geodesies in harmonic gauge. Those exact differential 
equations for the null geodesies will be solved numerically with high accuracy (see below) 
and that numerical solution provides the required reference. 



A. 



Metric tensor 



As it has been already discussed in Klioner & Zschocke [l| in harmonic gauge 








(1) 



dxf^ 
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the components of the covariant metric tensor of the Schwarzschild solution are given by 

1 — a 



fl'oo — ~ :r- — 5 
1 + a 

goi = 0, 



1 + a 



gij = (1 + a)^ % + ^ X* (2) 
x"' 1 — a 

where 

a = — , (3) 

a; 

m = ^ is the Schwarzschild radius of a body with mass M. The contravariant components 
read 

1 — a 

Considering that the determinant of the metric can be computed as 

9 = -{l + a)\ (5) 
one can easily check that this metric satisfies the harmonic conditions ([T]). 

B. ChristofFel symbols 



The Christoffel symbols of second kind are defined as 

2 ^ \dxP dx^ dx- ) ■ 

Using ([2]) and (jlj) one gets 

_a 1 , 
x2 l-a2 ^ ' 
, _ g 1 - « ^» 
(l + a)3'^ ' 

and all other Christoffel symbols vanish. 

C. Isotropic condition 

As it has been pointed out in Section 11. C of Klioner & Zschocke jl[ the condition of 
isotropy 



6 



leads to the following integral of the equations of light propagation 

where /z* is the coordinate direction of propagation (^t-^ = 1), as is the position of the photon 
and s is the absolute value of the coordinate light velocity normalized by c: s = |i|/c. 



D. Equation of isotropic geodesies 

Reparametrizing the geodetic equations 

by coordinate time t = (see e.g. Section 11. D of Klioner & Zschocke Q]) and using 
the Christoffel symbols computed above one gets the differential equations for the light 
propagation in metric ([2]): 



X - 



x^ 



a X ■ X 



—c 



[l + a) 



— X ■ X + a- 



1 -a2 



x + 2—- -{x-x)x. 11 

1 — a-^ 



Eq. ([9]) for the isotropic condition together with x ■ x = could be used to avoid the 
term containing i ■ i, but it does not simplify the equations. 



III. NUMERICAL INTEGRATION OF THE EQUATIONS OF LIGHT 

PROPAGATION 

Our goal is to integrate Eq. flTT]) numerically to get a solution for the trajectory of a light 
ray with an accuracy much higher than the goal accuracy of l/xas ~ 4.8 x 10~^^. For this 
numerical integration a simple FORTRAN 95 code using quadrupole (128 bit) arithmetic 
has been written. Numerical integrator ODEX [2| has been adapted to the quadrupole 
precision. ODEX is an extrapolation algorithm based on the explicit midpoint rule. It has 
automatic order selection, local accuracy control and dense output. Using forth and back 
integration to estimate the accuracy, each numerical integration is automatically checked 
to achieve a numerical accuracy of at least 10~^^ in the components of both position and 
velocity of the photon at each moment of time. 

The numerical integration is first used to solve the initial value problem for differential 
equations (ITTi) . Eq. ^ should be used to choose the initial conditions. The problem of 
light propagation has thus only 5 degrees of freedom: 3 degrees of freedom correspond to 
the position of the photon and two other degrees of freedom correspond to the unit direction 
of light propagation. The absolute value of the coordinate light velocity can be computed 
from ([9]). Fixing initial position of the photon a;(to) and initial direction of propagation /z 
one gets the initial velocity of the photon as function of n and s computed for given /it and 
x: 



x(to) = Xo, 

x{to) = cs 



(12) 
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The numerical integration yields the position x and velocity a; of a photon as function of 
time t. The dense output of ODEX allows one to obtain the position and velocity of the 
photon on a selected grid of moments of time. Eq. ([H]) holds for any moment of time as 
soon as it is satisfied by the initial conditions. Therefore, ([9]) can be also used to estimate 
the accuracy of numerical integration at each moment of integration. 

For the purposes of this work we need to have an accurate solution of two- value boundary 
problem. That is, a solution of Eq. ffTTl) with boundary conditions 

x{to)=XQ, 

x{t) = x, (13) 

where Xq and x are two given constants, to is assumed to be fixed and t is unknown and 
should be determined by solving flTTl) . Instead of using some numerical methods to solve 
this boundary problem directly, we generate solutions of a family of boundary problems 
from our solution of initial value problem (fT2l) . Each intermediate result computed by 
ODEX during the integration with initial conditions (|T2i) gives us a high-accuracy solution 
of the corresponding two- value boundary problem (|T3|) : t and x are just taken from the 
intermediate steps of our numerical integration. 

In the following discussion we will compare predictions of various analytical models for 
the unit direction of light propagation n{t) for a given moment of time t. The reference 
value for these comparisons can be derived directly from the numerical integration as 

n{t) = (14) 

The accuracy of this numerically computed n in our numerical integrations is guaranteed 
to be of the order of 10"^'* radiant and can be considered as exact for our purposes. 



IV. STANDARD POST-NEWTONIAN APPROACH 

In this Section we will recall the standard post-Newtonian approach and will compare the 
results for the light deflection with the accurate numerical solution of the geodetic equations 
described in the previous Section. 



A. Equations of post-Newtonian approach 

The well-known equations of light propagation in first post-Newtonian approximation 
with PPN parameters have been discussed by many authors. The differential equations for 
the light rays are given by the post-Newtonian terms of Eq. (22) of Klioner & Zschocke [l[: 



ax , , ax [x^ x^^ 



X 



(c^ + 7 x'^ x'^) ^ + 2 (1 + 7) ^ + 0(c-2) . (15) 

The analytical solution of f[T^ can be written in the form 

a^(t)=a;pN + 0(c-^), (16) 

iCpN = + c (t - to) cr + Aa;(i:) , (17) 
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where 

X + (T ■ X 



Ax(t) = — (1 + 7)m cr x (xq x cr) + cr log 

\ \x — cr ■ X Xq — cr • Xq J Xq + cr ■ Xq 

(18) 

Solution ffTBl) - ffT51) satisfies the following initial conditions: 

X{tQ)=XQ , 

lim x{t) = ca. (19) 

t— oo 

From Eqs. (fTB]) - ffT5]) it is easy to derive the following expression for the unit tangent vector 
at observer's position (note, in boundary problem we consider x^^ as the exact position cc, 
according to Eq. flTB]) ): 

, , . X [Xr\ X X) . . 

npN = fc- 1 + 7 ^, 20 

X\XXq\X-Xi^) 

where R = x — xq, k = R/R. By means of Eq. fj47|) given below we obtain that for the 
angle (5(npN, k) between ripN and k one has (for 7 = 1) 

where 

d = k X {xq X k) = k X {x X k) . (22) 
In the limit of a source at infinity one gets 

4:771 

lim 5(npN,fc)< — . (23) 



B. Comparison between the post-Newtonian approximation and numerical solution 



In order to determine the accuracy of the standard post-Newtonian approach we have 
to compare the post-Newtonian predictions of the light deflection with the results of the 
numerical solution of geodetic equations. Here, we compare the difference between the 
unit tangent vector defined by (|2U]) and the vector n calculated from the numerical 
integration using (JH]). 

Having performed extensive tests, we have found that, in the real solar system, the error 
of ripN for observations made by an observer situated in the vicinity of the Earth attains 16 
/xas. These results are illustrated by Table [T] and Fig. [H Table [T] contains the parameters 
we have used in our numerical simulations as well as the maximal deviation between ripN 
and n in each set of simulations. We have performed simulations with different bodies of 
the solar systems, assuming that the minimal impact distance d is equal to the radius of 
the corresponding body, and the maximal distance x between the gravitation body and the 
observer is given by the maximal distance between the gravitational body and the Earth. 
The simulation shows that the error of ripN is generally increasing for larger x and decreasing 
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for larger d. The dependence of the error of ripN for fixed d and x and increasing distance 
between the gravitating body and the source at Xq is given on Fig. [T] for the case of Jupiter, 
d being taken to be minimal and x to be maximal as given in Table [B Moreover, the error 
of npN is found to be proportional to which leads us to the necessity to deal with the 
post-post-Newtonian approximation for the light propagation. 























/' 








/ 









50 100 150 200 

Xq [AE] 



FIG. 1: The angle between ripN and n for Jupiter. The vector npN is evaluated by means of the 
standard Newtonian formula ()20p . while n is taken from the numerical integration as described in 
Section IIIII Impact parameter d is taken to be the radius of Jupiter and the distance x between 
Jupiter and the observer is 6 AU. 





Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


m = GM/c^ H 
dmin [10^ m] 

3^ max [AU] 


1476.6 
696.0 
1 


1476.6 
105781.7 
1 


1.40987 
71.492 
6 


0.42215 
60.268 
11 


0.064473 
25.559 
21 


0.076067 
24.764 
31 


<^max [/ias] 


3187.8 


6.32 X 10"^ 


16.13 


4.42 


2.58 


5.84 



TABLE I: Numerical parameters of the Sun and giant planets are taken from [J, |5|. dmm is the 
minimal value of the impact parameter d that was used in the simulations. For each body dmin 
are equal its radius. For the Sun at 45° the impact parameter is computed as d = sin 45° x 1 AU. 
^^max is the maximal absolute value of the position of observer x that was used in the simulations. 
Smax is the maximal angle between npN and n found in the numerical tests. 



V. POST-POST-NEWTONIAN SOLUTION OF BOUNDARY PROBLEM 

In [l| an explicit analytical solution of the parametrized post-post-Newtonian equations 
of light propagation in the gravitational field of one spherically symmetric static body has 
been derived. The solution a^ppN is given by Eqs. (26)-(35) of [l|. Boundary problem f|T3|l 
has been considered in In this Section, we derive analytical upper estimates of all the 
terms in the post-post-Newtonian solution of the boundary problem and find which terms 
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are responsible for numerical errors of the post-Newtonian solution described in the previous 
Section. In this way we derive the simplest possible formulas that agree with exact solution 
at a given numerical level. 



A. Analytical estimates of the individual terms in cr 



The propagation time between Xq and x is given by Eq. (50) of p|: 



N 



pN 



ApN 



ppN 



ppN 



CT=R 



+ (1 + 7) m log 



+ ^ (1 + 7)2m2 



X + Xq + R 
X + Xq — R 

R 



\X X Xq\ 



((x-xo)'-i?') 



1 

+ - ae „ 
8 R 



■m? f Xq — x'^ — R^ — Xq — R^ 



+ 



x^ 



Xn 



+ \ a (8(1 + 7) - 4/? + 3e) , ^ , 6{x, x^) 
4 \x X. xq\ 

+0(c-6) . 



(24) 

Here and below we classify the nature of the individual terms by labels N (Newtonian), pN 
(post-Newtonian), ppN (post-post-Newtonian) and ApN (terms that are formally of post- 
post-Newtonian order, but may numarically become significantly larger than other post- 
post-Newtonian terms, see below). Using \x x Xq\ = Rd where d is the impact parameter 
defined by (12^ . and assuming general- relativistic values of all parameters a = /5 = 7 = e = l 
one gets the following estimates of the sums of the terms labelled "ppN" and "ApN" : 



c5t = |c5r, 



ppN 



15 m 

<—TX — 

-4 d 



(25) 
(26) 



Estimates fl25|) - (l26|) are proved in Appendix |Al Note that here are below all estimates we 
give are reachable for some values of parameters and, in this sense, cannot be improved. 
From these estimates we can conclude that among the post-post-Newtonian terms c^tapn 
can become significantly larger compared to the other post-post-Newtonian terms. 

A series of additional Monte-Carlo tests using randomly chosen boundary conditions have 
been performed toverify the given estimates of the post-post-Newtonian terms. The results 
of these simulations are described in Table IIIII 

The effect of c6t for the Sun is less than 3.8 cm for arbitrary boundary conditions. 
Therefore, the formula for the time of light propagation between two given points can be 
simplified by taking only the relevant term: 

c r = it + (1 + 7) m log 

X + Xq — R 

+ 7)^ {R' -{X- XQ^) + (^) + 0{m% (27) 
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This expression can be written in an elegant form 



CT = R + (1 + 7) m loj 



X + xo + -R + (1 + 7) m 
X + xo — -R + (1 + 7) m 



+ 



m 



(28) 



that has been aheady derived by [3(] in an inconsistent way (see Section 8.3.1.1 and Eq. (8- 
54) of [3I). As a criterion if the additional post-post-Newtonian term is required for a given 
situation, one can use Eq. (l25l) giving the upper boundary of the additional term. 



B. Analytical estimates of the individual terms in transformation from fc to <t 

Transformation between k and cr is given by Eq. (51) of [l|: 

(T =k 



pN 



ApN 



scaling 



ppN 



ppN 



ppN 



-f (1 7) m — k X (xqX x) 

\x X aJoT 

1 1 
H — (1 + 7)^ A; X (xq x x) 



2 " ' ~ ' \x X xq\ 

:i + 7)' 2(x-Xo + i?)2 



(x -|- xo) (x — xo — R) (x — xo + Ry 



m 



\x X a^oP 



k 



+ 171^ k X (xo X x) 



1 1 



1 



xf 



+ - (8(l-h7-a7)(l + 7) -4a/5-F3ae) 



X X Xo\ 



2i?2 (tt - 6{k, x)) + (x^ -xl- R^) 5{x, x^)^ 



+0{c 



-6\ 



(29) 



Let us estimate the magnitude of the individual terms in Eq. 0291) in the angle 5(cr, k) 
between cr and k. This angle can be computed from vector product k x cr, and, therefore, 
the term in fl29p proportional to k and labelled as "scaling" plays no role. Here and below 
terms proportional to k do not influence the directions in the given order of magnitude, but 
are only necessary to keep the involved vectors to have unit length. The total effects of the 
terms of the other groups on x cr can be estimated using 



k X [k X {xo X x)]\ = \k X {xq x x) 



Rd. 



(30) 
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and general- relativistic values of the parameters a = /3 = 7 = e = l (see Appendix [B|) : 



IPpN 



< 



4m 



Xo < X, 



'' 4m 
X < —r- 

Xq > X a 



_ X + Xo 



IPapnI <16 



m 



27 



(x + Xo), I X < Xo < X, 



X^ Xo 
(X + Xof 



Xo < 9 X or Xo > X . 



P = IPppnI < 



15 m^ 



(31) 

(32) 
(33) 



Note that and p^pN themselves as well as their estimates are not continuous for x ^ Xq 
since in this limit an infinitely small change of x leads to big changes in k. Discontinuity 
of the same origin appears for many other terms. The limit x ^ Xq and the corresponding 
discontinuity have, clearly, no physical importance. 

We see that among terms of order m^ only IpapnI cannot be estimated as const x iv? /(P. 
The sum of the three other terms can be estimated as given by (1331) . The values of p for 
solar system bodies are given in Table [Tll In most cases these terms can be neglected at 
the level of 1 /las. Indeed, it is easy to see that p can be comparable with 1 /ias and even 
exceed this limit only for observations within 5 angular radii from the Sun. Again, Monte- 
Carlo simulations have been performed to check the actual maximal magnitude of these 
terms. The results are given in Table IIIII Accordingly, we obtain a simplified formula for 
the transformation from fc to cr keeping only the post-post-Newtonian term that can become 
larger than 1 pas also far from the Sun: 



fc + (1 + 

, (1 + 7)^ 



X - Xo + -R , , . 
7j m — — k X (a^o x x) 



\X X Xq\ 



m^ (x + Xo) 



(x — Xo + Ry (x — Xo — R) 

\x X Xf)\'^ 



kx{xoxx) + 0(^] + 0{m^' 



m 



This can be also written as 



CT = k +dS {l-S^{x + xo) + 



c n , ^ / 1 a;o -X 



m 



(34) 

(35) 
(36) 



where d is defined by f l22]) . Eq. (1521) can be used as a criterion if the additional post-post- 
Newtonian term in fl34p or (1551) is necessary for a given accuracy and configuration. 
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C. Analytical estimates of individual terms in transformation from <t to n 

Transformation between n and cr is given by Eq. (55) of [l[: 



N 



pN 



n = cr 



[1 + 'j)mk X {xq X x) 



R 



X X Xq 



1 + 



k ■ X 



X 



scaling 



ApN 



ppN 



ppN 



ppN 



ppN 



+ ^(l + 7)2m2 
k X (xq X X 



k 



+ (1 + 7)' 



X X Xq\^ X 

(1 + 7)' 
R 



^ 1 + ^— ^ 1 (3a; -xq- R){x-xo + R) 



X 



R 



X X Xq 



1 



X X Xo\ 

k ■ x\ 1 



1 + 



k ■ x\ R {R^ — (x — xo) 



2 I x X ajn P 



X X 



1 k ■ X 

-- ae 7 

2 Rx^ 



R 



\ (8(l + 7-a7)(l + 7) -4«/3+3ae)^|^^^^|, 



+0(c 



^1 + 7 - a7)(l + 7) -4a/5 +3ae) 



a; X a:;o 



(tt — 5(fc, a;)) 



(37) 



Let us estimate the magnitude of the individual terms in Eq. (1371) in the angle 5(cr, n) 
between n and cr. This angle can be computed from vector product cr x n, and, therefore, 
the term in (1371) proportional to k plays no role since a xk = Oim). To estimate the effects 
of the other terms in (1371) we take into account that 



(7 X [kx (xo X a;)) | = Rd^ 0{ 



m) 



and assume again a = /? = 7 = e = l. We get (see Appendix ICj) 



\cpp^\=2m 



I ApN I 



cr X [fc X {xo X X 
(7 X [k X {xq X x)] 



R 






X X Xo 


2 





m 



(38) 



(39) 



1 



k ■ X 



X 



- (x - Xof 



X X Xq 



< 4 , 

d"^ d {x + xo)^ (i^ d ' 



. 15 
<^ = |¥'ppnI ^ rf^- 



(40) 
(41) 



Again here y'ppAr is the sum of all the terms in (157|) labelled as "ppN" . These terms can attain 
1 nas only if one observes within 5 angular radii from the Sun. The values of ip for solar 
system bodies are given in Table [Tll In most cases these terms can be neglected at the level 
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of 1 /ias. Again, Monte-Carlo simulations have been performed to check the actual maximal 
magnitude of these terms. The results can be found in Table IIIII Accordingly, we obtain a 
simplified formula for the transformation from cr to n keeping only the post-post-Newtonian 
term that can become significantly larger than the others: 



N d f ^ k ■ X 

n = (T - (l + 7)m— IH 

\ X 



d 



k ■ x\ B? — [x ~ xq) 



2R 



m 



+ 0{'—]+ 0{m^' 



where d is defined by (1221) . The same formula can be written as 



R + Xo-x \ ^ ^ /m^ 
R + Xq + X J yd"^ 



n = (T + dT I 1 + Tx 

T = 



C(m^) 



, m / k ■ X 



(42) 



(43) 
(44) 



D. Analytical estimates of individual terms in transformation from to n 

Transformation between n and cr is given by Eq. (56)-(57) of Q: 



N 



pN 



ApN 



scaling 



ppN 



ppN 



ppN 



ppN 



where 



n = k 



-(1 + 7) m ■ 



k X {xq X x) 

X {x Xq + X ■ Xq] 



kx (xoxx) 
-(1 + 7) m — F 

X [X Xq + X ■ Xo) 
1 m2 ((x - xo)2 - 



--(l + 7)^^fc 



x^ 



\X X Xq\ 



m k X [xq X x) 



1 1/1 
+ - ae — 



(1 + 7)' 



R^ — {x — Xq) 

\x X Xq\^ 



4:'~~R\Rxl Rx^ 



k ■ X 



x" 



-\ (8(l + 7-a7)(l + 7) -4a/? + 3ae) R ^ ' ^ ,„ 
4 x^\xx xq\^ 

+ \ (8(l + 7-a7)(l + 7)-4a/3 + 3ae) '"^'^ ~ f 5{x,Xo) 



\X X Xq\ 



+ O (c-6) 



F = -(1 + 7)m 



X + Xq 



XXq + X ■ Xq 



(45) 



(46) 
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As in other cases our goal is to estimate the effect of the individual terms in Eq. fl45p on the 
angle 6{k, n) between k and n. This angle can be computed from vector product k x n. 
The term in (H^ proportional to k obviously plays no role here and can be ignored. For 
the other terms taking into account Eq. fl30|) and considering the general-relativistic values 
a = /3 = 7 = e = l one gets (see Appendix [D]) 

1 Ifc X (xn X x)\ m Xn m , . 

o^pN =2m- J ^ < 4-^^ < 4-, 47 

X xxq + X ■ xq a X + xq a 

X X Xq + X ■ Xq 

^„m? RxXn -,„rn'^ xxn _^m^x 



or, alternatively, 



, 64 



We give four possible estimates of |<-l»apn|- These estimates can be useful in different situa- 
tions. Note that the last estimate in fHHl) and the estimate in (1491) cannot be related to each 
other and reflect different properties of Iu^apnI as function of multiple variables. 

The sum of all the terms in (H^ labelled as "ppN" is denoted as uJppN and can be 
estimated as 

15 

^ = |^ppn| < — TT — . (50) 

Again these terms can attain 1 //as only for observations within 5 angular radii from the Sun. 
The values of u for solar system bodies are given in Table [Tll One can see that these terms 
can be neglected at the level of 1 /las in most cases. The corresponding results of our Monte- 
Carlo simulations can be found in Table IIIII Accordingly, we obtain a simplified formula 
for the transformation from to n keeping only the terms which cannot be estimated as 

n = fe-(l + 7)ml ^ (1 + F) + O f ^ V ^i^') ' (51) 

X xxo + X ■ xq \d J 

where F is given by fH^ . This can be also written as 

Xq + x\ ^ . z\ 



n = fc + dP + +0 J +0(m^), (52) 

^ m / Xq — X k ■ x\ , . 

where d is given by Eq. fl22l) . Let us also note that the post-post-Newtonian term in flSTj) 
and fl52|) is maximal for sources at infinity: 



k^ApN < lim cjApN = lim 1 + 7 m-^ ^— '-^ F = {1 + -ff {1 - cos<^y — - , 

xo^oo xo^oo X XXq + X-Xq d'^ a 

(54) 
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c" 0.02 

^ a 
c 



.00 

50 100 150 200 

Xq [AE] 

FIG. 2: The angle between n'^^ and n for Jupiter. Vector n'^^ is evaluated with the aid of 
()5ip . while n is the high-accuracy numerical solution of exact geodetic equations as described in 
Section IIIII Impact parameter d is taken to be the radius of Jupiter and the distance x between 
Jupiter and the observer is 6 AU as on Fig. [TJ This figure demonstrates the numerical validity of 
(j52p at the level well below 1 /ias. 





Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


c6t [10-^ m] 


36906.0 


242.9 


0.328 


0.036 


0.002 


0.003 


p, cp, Tp, LO [10^^ fias] 


10937.4 


0.474 


0.945 


0.120 


0.016 


0.023 



TABLE II: Numerical values of the analytical upper estimates of the post-post-Newtonian terms of 
order oiO{m^/d) in Eq. ([MD of order 0{m'^/d'^) in Eqs. 1^, ([37]), (gS]), and dSS]). The analytical 
estimates are given by Eqs. ([26]) . (f33|) . (jH]), ([50]) . and ([6T|) . respectively. One can see that at the 
level of 10 cm in the distance and 1 fias in angles these terms are irrelevant except for observations 
within 5 angular radii from the Sun. 



where $ = S{xo, x) is the angle between vectors Xq and x. Several useful estimates of these 
terms are given by fH8l) - (l49l) . These estimates can be used as a criterion which allows one 
to decide if the post-post-Newtonian correction is important for a particular situation. 

Using estimate (15^ and the parameters of the solar system bodies given in Table [T] one 
gets the maximal values of the post-post-Newtonian correction shown in Table IIVI For 
grazing rays one can apply cos*!* ~ —1, while for the Sun at 45° one can apply cos$ ~ 
— 1 / -\/2- Comparing these values with those in the last line of Table [T] one sees that the 
post-post-Newtonian correction matches the error of the standard post-Newtonian formula. 
The deviation for a grazing ray to the Sun is a few /xas and originates from the post-post- 
Newtonian terms neglected in Eq. fl52|) . Vector n computed by fl52l) can be denoted as 
(a post-Newtonian formula enhanced by one post-post-Newtonian term that can become 
large). The numerical validity of n'^j^ can be confirmed by direct comparisons of n'^j^ and 
vector n computed by numerical integration of geodetic equations as discussed in Section 
IIVBI The results for Jupiter are given on Fig. [2] (cf. Fig. [T]). 
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Sun 


Qiin of 

oun at 40 


Jupiter 


Saturn 


Uranus 


Neptune 


c8t [10-6 m] 


36846.9 


178.71 


0.327 


0.0348 


0.00192 


0.00274 


p [10^^ X //as] 


10747.6 


0.349 


0.942 


0.119 


0.0154 


0.0228 


if [IQ^^ X /ias] 


10713.2 


0.112 


0.942 


0.119 


0.0154 


0.0228 


^ [10^^ X ^as] 


10713.2 


0.112 


0.942 


0.119 


0.0154 


0.0228 


UJ [lO^-^ X ^as] 


10682.8 


0.239 


0.834 


0.096 


0.0109 


0.0134 



TABLE III: Maximal values of the sum of the terms of order of 0(m^/d) in Eq. ()24p and of order 
0{rn?' /(f') in Eqs. ()29p . ()37p . ()45p and (|58p obtained from numerical simulations. Two simulations 
have been performed. For the first simulation 10^ starting points xq were taken within 50 AU from 
the relevant massive body. For the second simulation 10^ starting points were taken at a random 
distance but with a constraint that in each case the straight line between the starting and final 
points is tangent to the surface of the body under consideration. Final point x is always chosen on 
the orbit of the Earth. The position on the Earth orbit is taken randomly. For each of the 2 x 10^ 
points the corresponding terms were evaluated numerically and the maximal value is given in the 
Table. The fact that these values are always below the corresponding analytical estimations given 
in Table [TI] can be considered as additional confirmation of the estimates ([26]) , (j33]) , ([^T]) , ([50]) , and 
(|6ip . Note that the values of oj for all cases and all the values for "Sun at 45°" are systematically 
smaller than the estimates given in Table [III This behaviour is well understood and expected in 
the described set-up of the Monte-Carlo simulations. 





Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


max o^ApN [/^as] 


3192.8 


0.663 X 10-3 


16.11 


4.42 


2.58 


5.83 



TABLE IV: Maximal numerical values (|54p of the post-post-Newtonian correction in Eqs. ()5ip 
or ()52p for the solar system bodies given in Table [H A comparison of these values with the last 
column of Table |T] allows one to conclude that the post-post-Newtonian term in (j5ip is responsible 
for the errors of the the standard post-Newtonian formula ([20]) . 



VI. TRANSFORMATION FROM fc TO n FOR STARS AND QUASARS 

In principle, the formulas for the boundary problem given above are valid also for stars 
and quasars. However, for sufficiently large xq the formulas could be simplified. It is the 
purpose of this Section to derive necessary formulas for this case. 

A. Transformation from k to cr 

First, let us show that for stars and quasars the approximation 

cT = k (55) 
is valid for an accuracy of 1 /las. Using estimates (131]) and (l32j) for the two terms in Eq. 
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Xf) [pel 


Sun 


Sun at 45° 


Jupiter 


Saturn 


Uranus 


Neptune 


1 


8.506 


0.056 


0.473 


0.309 


0.212 


0.382 


10 


0.851 


5.586 X IQ-'^ 


0.047 


0.031 


0.021 


0.038 


100 


0.085 


0.559 X 10-3 


4.740 X 10-3 


3.086 X 10-3 


2.122 X 10-3 


3.819 X 10-3 



TABLE V: Numerical values of estimate ()56p in /ias for the light deflection due to the solar system 
bodies with various values of xq. 

( 134|) one can see that for xq ^ x the angle 6{(t, k) can be estimated as 

S((T k) < 4^ ^ ( 1 + 4^ ^ ^° ^ (56) 

' ~ d X + Xq \ d d X + Xq J 

Clearly, S{(T,k) goes to zero for xq —>■ oo. Numerical values of this estimate are given in 
Table IVl for xq equal to 1, 10 and 100 pc. Angle S{cr, k) is smaller for stars at larger distances. 
However, for hypothetical objects with xq < 1 pc the difference between cr and k must be 
explicitly taken into account. 

B. Transformation from cr to n 

As soon as we accept the equality of cr and k for our case the only relevant step is 
the transformation between cr and n. This transformation in the post-post-Newtonian 
approximation is given by Eqs. (53)-(54) of Introducing impact vector computed using 
cr and the position of the observer x 



= cr X (x X a) 



(57) 
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we can re-write Eqs. (53)-(54) of fi\ as 



N 



pN 



ApN 



scaling 



ppN 



ppN 



ppN 



ppN 



n 



+ (1 + 7)^m 



2 _2"'o' 



(T ■ X 

f 

X 



1 . 
-- m'a e 



a ■ X 



x" 



+ (1 + 7)'^'^- (1 + 
di X 



X 



<j ■ X 



X 



(8(l+7-a7)(l + 7)-4«/3 + 3«e) ^ ^-J^ 

(8(1 + 7 -a7) (1 + 7) -4a/5 + 3ae) (tt - 5(<t, a;)) 

4 d^ 

+0(m='), 



(58) 



where dn 



\dn 



\<7 ^ x\. Now we need to estimate the effect of the individual terms in 
Eq. 0581) on the angle (5(cr, n) between a and n. This angle can be computed from vector 
product cr X n. The term in (1581) proportional to a obviously plays no role and can be 
ignored. For the other terms taking into account that |cr x cZo-| = do- and considering the 
general-relativistic values a = /3 = 7 = e = lwe get 

Icr X d„ 



pNl 



2m- 



dl 



_^ (7 ■ X\ ^ ^ 

X J ~ d„' 



1^ 



ApN I 



Am 



Icr X d„| X 



dl 



dn 



cr ■ X 



X 



m^ X 



< 16 

dl d„ 



I / I 15 
^ = IV'ppnI < d^' 



(59) 
(60) 
(61) 



where V'ppN is the sum of all terms of order m^/d^. in (1581) . Estimate (1611) obviously agrees 
with estimate (1411) for ip. Numerical values of this estimate can be found in Table [ITl 
The estimates show that these terms can be neglected at the level of 1 yuas except for 
the observations within 5 angular radii from the Sun. Omitting these terms one gets an 
expression valid at the level of 1 /las in all other cases: 

(62) 



m 



n = (T + d„Q{l + Qx)+0^j^ ) +C(m3), 
^ , ^ m / ax 



(63) 



Note that for xq ^ oo this coincides with f lS2]) - flS5]) and with fH51) - fH^ . This formula 
together with cr = k can be applied for sources at distances larger than 1 pc to attain the 
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accuracy of 1 //as. Alternatively Eqs. can be used for the same purpose giving 

slightly better accuracy for very close stars. However, distance information (parallax) is 
necessary to use (132]) - (155]) . 

VII. SUMMARY AND CONCLUDING REMARKS 

In this report the numerical accuracy of the post-Newtonian and post-post-Newtonian 
formulas for light propagation in the parametrized Schwarzschild field has been investi- 
gated. Analytical formulas have been compared with high-accuracy numerical integrations 
of the geodetic equations. In this way we demonstrate that the error of the standard post- 
Newtonian formulas for the boundary problem (light propagation between two given points) 
cannot be used at the accuracy level of 1 yuas for observations performed by an observer 
situated within the solar system. The error of the standard formula may attain ~ 16 /ias. 
Detailed analysis has shown that the error is of post-post-Newtonian order 0{m?). On the 
other hand, the post-post-Newtonian terms are often thought to be of order /(P and can 
be estimated to be much smaller than 1 /xas in this case. To clarify this contradiction we 
have investigated the post-post-Newtonian solution for the hght propagation derived in 
For each individual term in relevant formulas upper estimates have been found. It turns out 
that in each case one post-post-Newtonian term may become much larger than the other 
ones and cannot be estimates as const x iv? /d?. These terms depend only on 7 and do not 
come from the post-post-Newtonian terms of the corresponding metric. The formulas for 
transformations between directions cr, n and k containing both post-Newtonian terms and 
post-post-Newtonian ones that can be relevant at the level of 10 cm for the Shapiro delay 
and 1 /ias for the directions have been derived. The formulas are given by Eqs. f l28p . fl35p - 
dSS]), (llSD-dllD, dSl-dM]), and (in2)-(inSD- These formulas should be considered as formulas 
that guarantee this numerical accuracy. 

The derived analytical solution shows that no "native" post-post-Newtonian terms are 
relevant for the accuracy of 1 yuas in the conditions of this note (no observations closer than 
five angular radii of the Sun). "Native" refers here to the terms coming from the post-post- 
Newtonian terms in the metric tensor. It is, therefore, not the post-Newtonian solution 
itself, but the standard analytical way to convert the solution of the initial value problem 
into the solution for the boundary problem that is responsible for the numerical error of 16 
/xas mentioned above. 

Let us finally note that the post-post-Newtonian term in fl321) - flS^ is closely related to 
the standard gravitation lens formula. Here we only note that all the formulas given in [l| 
and in this paper are not valid for d = Q {d always appear in the denominators of the relevant 
formulas). On the other hand, the standard post-Newtonian lens equation successfully treats 
this case, known as the Einstein ring solution. The relation between the lens approximation 
and the standard post-Newtonian expansion is a different topic which will be considered in 
a subsequent paper. 
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APPENDIX A: ESTIMATES OF TERMS IN THE SHAPIRO DELAY 

In order to get fl2^ we write the corresponding term as 



Ic^tapnI 



2 m' 



R 



\X X Xo\ 



{{x-xoy-R^) 



m 



2z{l - cos$) 
d? 1 + - 2z cos$ 



2 — R 



(Al) 



where $ = 5{x, Xq) is the angle between x and ajg, and z = Xq/x. It is easy to see that for 
0<<l><7randz>0 



_ 2z{l - cos<!>) ^ Az ^ ^ 
•^^ ~ 1 + ^2 _ 2^ cos $ - (1 + zY - 



(A2) 



This immediately gives fl25|) . Here and below we always give estimates that cannot be 
improved in the sense that they are reachable for certain values of the parameters. 
For fl26l) we write 



c6t = \c6t. 



PpNl 



1 /xg — — -R^ 



8 R 
1 

4~d 



X 



xl-R^ 



x^ 



Xn 



R 



15 2 

— m 

4 |a; X ajQ 



sin $ 



2;^ cos $ — 2 2; + cos $ 
1 + - 2 z cos $ 



15$ 



(A3) 



Here and below $ = 5(a;,a[;o) is the angle between x and Xq, and z = xq/x. One can 
show that for < $ < tt and z > 



f2 



z^ COS $ — 22; + cos $ 

sin $ h 15 $ 

1 + _ 2 z cos $ 



< 15 vr. 



(A4) 



and this immediately gives fl2Bl) . 



APPENDIX B: ESTIMATES OF TERMS IN THE TRANSFORMATION 

BETWEEN (7 AND k 



For Eq. (1311) we note that 



. . X — Xq + R ,^ , , . 2m X — Xo + -R 

IPpnI =2m -p— — P5- |fc X (ajo x a^J 



la; X Xq\ 



d 



R 



2m 



1-2 



d + z"^ - 22: COS $ 



+ 1 



(Bl) 



Again $ = 5(a;, a^o) is the angle between x and xq, and 2; = xq/x. One can show that for 
0<$<7rand;z>0 



fs 



1-z 



Vl + -2^ - 22; COS $ 



1< 



Z < 1 

2 <2 

2; > 1 



(B2) 



L 1 + z' 



and this leads to fl^ . The discontinuity of and its estimate at z ^ 1 are discussed in 
the main text after Eq. 
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The term p^pN can be written as 



m 



|pApNl = 2^(a; + xo) 



{x-xo + R) [{x - xof - ] 



i?3 



,rn? X , , 1 — cos $ 
(P d ^ ' l + z^-2z cos$ 



1 + 



1 - z 



Vl + -2^ - 2^; cos $ 



(B3) 



For < $ < TT and 2; > one has 
l-cos$ 

/4 = ^ (1 + ^) 



1 + 



l-Z 



1 + ^2 _ 2z COS $ V Vl + ^2 _ 2^ cos$ 



< 



if(l + ^), i<^<l, 

z 

(B4) 

This gives Eq. (1321) . The function itself and its estimate flB4p are again not continuous for 
z = 1 (implying x = Xq). This is discussed after Eq. f l32|) . 
In order to get (1551) we write 



P = IPppnI ="^^ |fc X (^0 X x) 



1 1 
4 ^ 



1 1 



X^ XT 



15 



1 



X 2i?2 



(5(fc,aj)) + (, 



X 



1 

11? 



z{z^ - 1) sin^ $ 
(1 + z2 _ 2z cos$)^ 

z (cos $ - 2;) $ 



— 15 arccos 



8 |a; X a^oP 
xl - R^) 5{x,xq) 
1 — z cos $ 



2^ 



2 2 COS $ 



+ 15 



1 + 2;2 - 2 2; cos $ 



15 TT 



(B5) 



One can show that for < $ < vr and z > 

z{z^ - 1) sin^ $ 



(1 + ^2 - 2^ COS$)^ 



— 15 arccos 



1 — z cos $ 



2 (cos $ - 2;) $ 

1 + - 2 2; cos $ 



Vl + - 22; cos<l> 
< 157r. 



(B6) 



This immediately leads to fl55]) . 
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APPENDIX C: ESTIMATES OF TERMS IN THE TRANSFORMATION 

BETWEEN n AND a 



Estimate f l39D for cp^^ is trivial. For estimate ( HOl) of g?Ar,M we write 



^ApN 



(7 X [k X {xq X x)] 



1 + 



k ■ X 



^rn? R _ 

d? d \ X 



X / \X X Xq 

k ■ x\ R^ — {x — xq^ 



2i?2 



R 

4 — — 



l-zcos$ 2:(l-cos<l>) 



d"^ d \ + z"^ - 2z cos<l> y 1 + z2 - 2z cos $ 



(CI) 



where again $ = (5(a:;, a^o) is the angle between x and Xo, and z = xo/x. It is easy to see 
that for < $ < TT and z > 



1 — Z COS $ 



^ -cos'l') ^ 4z ^ ^ 



Vl + z^ -2z cos$ y l + z2-2zcos$ (1 + z)^ 
This immediately leads to fHU]) . For Eq. fHTl) we write 



(C2) 



= Iv'ppnI ="^^|fc X (^0 X 



15 k ■ X 



R 



X X Xq 

R 



k ■ x\ 1 1 k ■ X 



X / X 

15 R^ 



4 I a; X P 4 | a; x a^o 



2 

[tt — (5(fc, x)) 



1 

4 ^ 

1 
4 



16 -H 2 15 (vr - x)) 

XXX V X / X 



16 sin + sin cos - 2 sin"* \E' cos - 15 vr + 15 \I' 



(C3) 



where \& = S{k, x) is the angle between vectors k and x. Here we used that k ■ x = x cos \E' 
and d = Ifc X a;| = X sin \E'. For < \& < tt we have 



16 sin + sin vE^ cos \E' - 2 sin'^ cos \E' - 15 tt + 15 ^ 



< 157r 



(C4) 



and this proves Eq. (HTj) . 



APPENDIX D: ESTIMATES OF TERMS IN THE TRANSFORMATION 

BETWEEN n AND k 

In order to get fH7j) we write 

1 Ifc X (ajn X a;)| mxxn — aj-ajn fn 2;(1 — cos$) , 

cUpN =2m- J ^— - = 2 — — - — = 2 — ^ ^ ' =. Dl 

X xxo + ai-a^o d xR a v 1 + 2;2 - 2^; cos $ 
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Here again $ = S{x^ Xq) is the angle between x and a^o, and z = Xq/x. One can show that 
for < $ < vr and z > 



fs 



z(l-cos$) ^ 2z 



VI + ^2 - 2z cos $ 1 + 2: 
that immediately gives (147|) . To derive (l48l) and (l49i) we write 



(D2) 



I^ApnI = 2m 
= 4 



1 |A; X (ajo X a;)| 2 I ^ x (*o x a?) | a; + xq 

\r \ = 4 771 



X a;a;o + a;-a;o ■ ■ (xxo + aj-ajo) 

R z'^{l + z) (1 - cos$)^ 



X 



(P d {l + z'^-2z cos^f ' 
For < $ < TT and z > one gets 



(D3) 



/9 



< 



4z' 



z^{l + z){l -cos<l>)^ 
(1 + ^2 -2^ cos$)^ ~ {l + z 



(D4) 



< 1 and 



< 



This gives the first estimate in (HHj) . Trivial inequalities R < x+xq, (J^^ ^ jj^:^ ^ 27 

give the second and third estimates in fHHj) and estimate fj49l) . respectively. 
For fISOll we write 



a; = I^ppnI =777^ I fc X (a^o x x) 



^ R'-jx-xof ^ 1 



la; X ajoP 



4/? Vi^xg i?x2 



• X 



x'' 



• X 



15 X2 - x2 - i?2 



4 x^ I a; X a;o P 



+ ^ 



8 \x X xo\ 



6{x,Xo) 



1 



z {IQz - z cos $ - 15) sin ^ z{l - d,z^ + 2z^ cos $) sin^ $ 



+ 



1 + z"^ -2z cos $ 
15^ (cos $ - 2;) $ 



(l + 22_22 cos<l>)' 



1 + _ 2z cos $ 



(D5) 



For < $ < TT and z > Q one can demonstrate that 

/lO 



2; (16z - z cos $ - 15) sin $ ^ 2;(1 - S^;^ + 2z'^ cos $) sin^ $ 15^; (cos $ - 2;) $ 



1 + z"^ -2z cos $ 

< 15 7r 
and this leads to fl50|l . 



(1 + _2z cos^y 



1 + z"^ -2z cos $ 



(D6) 
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APPENDIX E: ESTIMATES OF TERMS IN THE TRANSFORMATION 
BETWEEN n AND a FOR STARS AND QUASARS 



Estimates (l59D-(l60l) are trivial. For fl6TD we write 



lax 



1 

111 

1 
4 



- 2 



2 

cr ■ X dcr cr ■ X 



1 / cr ■ x\ 15 cr ■ X 15n — 5{(T,x) 



+ — + 16 — - 15 TT + 15 6(cr, x) 



16 sin '^^ + cos sin - 2 sin^ cos - 15 vr + 15 



(El) 



where = 5{cr,x) is the angle between vectors cr and a?. Here we use d^j = \a ^ x\ = 
X sin\E'o- and cr ■ x = x cos^'o-- Therefore, for < < tt one can use estimate flC4p for /y 
to prove (1^ . 



